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Abstract. Solutions of initial-boundary value problems for systems of conser- 
vation laws depend on the underlying viscous mechanism, namely different vis- 
cosity operators lead to different limit solutions. Standard numerical schemes 
for approximating conservation laws do not take into account this fact and 
converge to solutions that are not necessarily physically relevant. We design 
numerical schemes that incorporate explicit information about the underlying 
viscosity mechanism and approximate the physically relevant solution. Nu- 
merical experiments illustrating the robust performance of these schemes are 
presented. 



1. Introduction 

Many problems in physics and engineering are modeled by systems of conserva- 
tion laws 

(1.1) U t + F(U) s = 0. 

Here, U:!lxI + 4 W n is the vector of unknowns and F : R m -> W m is the flux 
vector. The spatial domain is a set Q C K. The above equations are augmented with 
initial data. If fl is a bounded domain, the conservation laws are augmented with 
suitable boundary conditions. Examples of conservation laws include the shallow 
water equations of oceanography, the Euler equations of gas dynamics and the 
equations of MagnetoHydroDynamics (MHD). 

We assume that the system of conservation laws is strictly hyperbolic, i.e. the 
eigenvalues of the Jacobian matrix Fu are real and distinct. Also, we assume that 
all the eigenvalues of Fu are bounded away from 0: 

(1.2) Ax(U) < • • • < A fc (U) < -d < < d < A fe+1 (U) < • • • < A m (U) 
for some positive constant d > and some integer k < m. 



It is well known (see Dafermos [8 :j Chapter 6] ) that in general solutions of (1.1) 



form discontinuities (shock waves, contact discontinuities) in finite time even when 



the initial data are smooth. Hence, solutions of (1.1) are defined in the sense of 
distributions. 

In general, the distributional solution of a given Cauchy or initial-boundary value 
problem is not unique and hence various admissibility conditions have been intro- 
duced in the attempts at selecting a unique solution, see the book by Dafermos [SJ 
Chapters 4 and 8] for an extended discussion. These approaches often involve the 
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celebrated entropy condition, which can be formulated as follows: assume that sys- 



tem ( 1.1 ) admits an entropy-entropy flux pair, namely there exists a convex function 



S : R m -> K and a function Q : R m -s- R such that 

(1.3) Qu = 5uFu, 

where Su and Qu denote the gradients of the function S and Q, respectively. A 
distributional solution U satisfies the entropy admissibility condition if the following 
inequality holds in the sense of distributions: 

(1.4) S(U) t + Q(U) a < 0. 

Here, two remarks are in order: first, in general physical systems admit entropy- 



entropy flux pairs. Second, systems of conservation laws like (1.1 1 are derived by 
neglecting small scale effects like diffusion. Inclusion of these small effects in (1.1) 
results in the mixed hyperbolic-parabolic system: 

(1.5) U t e + F(U% = e(S(lT)U^. 

Here, e is a (small) viscosity parameter and B : R m — > R mxm is the viscosity matrix. 
For example, the Navier-Stokes equations are a viscous regularization of the Euler 
equations of gas dynamics. In physical systems, the entropy admissibility criterion 
is consistent with the zero small scale effects limit, namely one can show that, if the 



then the limit satisfies (1.4) 



solutions of the viscous approximation (1.5) converge in a strong enough topology, 



We now focus on the initial-boundary value problem obtained by coupling the 



system of conservation laws ( 1.1 ) with the Cauchy and Dirichlet data 

(1.6) U(x,Q) =U (x), xen^iXuoo) U(X u t) = U(i), teR+. 

The study of the initial-boundary value problem poses additional difficulties as 



compared to the study of the Cauchy problem: first, the problem (1.1 (-(1.6 1 is, 
in general, ill posed (i.e. it possesses no solutions) unless additional conditions are 
imposed on the data U. Possible admissibility criteria on U are discussed in Dubois 
and LeFloch [5]. 

Another additional difficulty one has to tackle when studying initial-boundary 



value problems is the following: consider the viscous approximation (1.5) coupled 
with the initial and boundary data 

(1.7) U e (a;,0) = U (.t), x e n = (A,, oo) U £ (A ; , t) = U;(i), t 6 R + . 



Assume that the initial-boundary value problem (1.5), ( |1.7[ ) is well posed (this is 



not always the case in the case when B is singular) and that for e — > + the solutions 
converge in a suitable topology to a limit U. In general, because of boundary layer 
phenomena, U may not satisfy the boundary condition U;(t) pointwise. Dubois 



and LeFloch [9] showed that, if the solutions of the viscous approximation (1.5) 



converge as e — > + to a solution of the initial-boundary problem (1.1 1, (1.6 1 in a 
sufficiently strong topology, then the following inequality holds: 

(1.8) Q(U(t)) - Q(U,(t)) - (5u(U,(i)), (F(U(f)) - F(U,(t)))) < 0. 

Here (•, •) denotes the standard scalar product in R m . 

A further difficulty in the study of initial-boundary value problems was pointed 
out in the works by Gisclon and Serre [TH H3]: they showed that the limit of 



the viscous approximation (1.5) depends on the underlying viscosity mechanism. 
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In other words, the limit of (1.5) in general changes if one changes the viscosity 
matrix B. 



As an example, we consider the linearized shallow water equations (2.13) with 



initial data (2.18) and boundary data (2.19). The system is a linear, strictly hy- 



perbolic, 2x2 system and is the simplest possible problem that can be considered 
in this context. We consider two different viscosity operators: an artificial uniform 



(Laplacian) viscosity (2.15) and the physical eddy viscosity (2.14). The resulting 
limit solutions are shown in the left of figure [T] As shown in the figure, there is 
a significant difference in solutions (near the boundary) corresponding to different 
viscosity operators. 

An extended discussion concerning the initial boundary value problem for sys- 
tems of conservation laws and its viscous approximation can be found in the books 
by Serre [25, 26 , Chapters 14 and 15], while we refer to the lecture notes by Serre [57] 
and to the rich bibliography therein for the theoretical treatment of the discrete 
approximation of viscous shock profiles. To conclude, we stress that analytically 
establishing the convergence e —¥ + for (1.5) is still an open problem in the gen- 



eral case, but results are available in more specific cases: in particular, Gisclon |12j 
showed local-in-time convergence in the case when B is invertible and, by extend- 
ing the analysis in Bianchini and Bressan [BJ, Ancona and Bianchini [3] proved 
global-in-time convergence in the case when B is the identity. 



1.1. Numerical schemes. Numerical schemes play a very important role in the 
study of system of conservation laws. Conservative finite difference (finite volume) 
methods are among the most popular discretization frameworks for ( |1.1[ ). See 
the book by LeVeque |23) for an extended discussion. Given the real numbers 
Xi < X r , we discretize the computational domain [Xi, X r ] by N + 1 equally spaced 
points Xj+1/2 = Xi + jAx with X±/2 = Xi and with mesh size Ax and we set 
Xj = x J- 1 /2+ a; j+i/2 ^ Time is discretized with a time step At 71 . The mesh size and 
time step are related by a standard CFL condition. 

The aim is to approximate cell averages U™ of the unknown U in the cell Cj = 



7-1/2) ) at time t n by the scheme 



(1.9) 



US' 



At™ 
Ax 



( F™ 

V 3+1 



/2 



' 3-1/2 



Here, F™ +1 y 2 = F(U",U™ +1 ) is the numerical flux. The numerical flux is obtained 



u " +1 



by solving (approximately) the Riemann problem for (1.1 ) with the states U™ and 

X\ are imposed by 



Following [23] , the Dirichlet boundary conditions at X 
setting in the ghost cell [s-i/2, #i/2] : 



(1.10) 



Un 



U,(t n ). 



However, standard numerical schemes may not converge to the physically rel- 
evant solution of the initial boundary value problem for a system of conservation 
laws. We illustrate this by again considering the linearized shallow water equa- 
tions (2.13) with initial data (2.18) and boundary data (2.19). The results with a 
standard Roe (Godunov) scheme for this linear system are presented in figure [lj 
right. The figure clearly shows that the Roe scheme converges to a solution that 
is different from the physically relevant solution of the system, realized as a limit 
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(a) Viscous profile 



(b) Roe scheme 



Figure 1. Left: The limit viscous profile for the linearized shallow 



water equations (2.13) with uniform (Laplacian) viscosity (2.151 
and eddy viscosity (2.14). Right: Approximate solutions with the 



Roe (Godunov) scheme for the linearized shallow water equations 
for the same data as the left figure. 



of the eddy viscosity approximation (2.14). In fact, the solution converges to the 



limit of the artificial uniform viscosity approximation (2.15) 



The problem with standard numerical schemes approximating the initial-boundary 



value problem ( 1.1 ) lies in the fact that they do not incorporate explicit information 



about the underlying viscous approximation (1.5 1. The implicit numerical viscosity 



added by such schemes may lead to the schemes converging to an incorrect solution. 
This situation presents analogies with the numerical approximation of non-classical 
shocks (see LeFloch[19 ), non-conservative hyperbolic systems (see Castro, LeFloch, 
Munoz Ruiz and Pares [7]) and conservation laws with discontinuous coefficients 
(see Admiurthi, Mishra and Veerappa Gowda [5]). 

Here, we design numerical schemes that incorporate explicit information about 
the underlying viscous operators. Consequently, these schemes approximate the 
physically relevant solutions of system of conservation laws. The schemes are based 
on the following two ingredients: 



(i.) An entropy conservative discretization of the flux F in ( 1.1 ) (see Fjordholm, 

Mishra and Tadmor [TT] and Tadmor [25]). 
(ii.) Numerical diffusion operators for ( |1.1[ ) that are based on the underlying 

viscosity matrix B in (1.5). 



We present both first- and second-order schemes that are shown (numerically) to 
converge to the physically relevant solution of the system of conservation laws. 

The rest of the paper is organized as follows: in Section |2j we discuss the the- 
oretical results concerning the initial-boundary value problems (1.5) that will be 



used in the following sections. In particular, explicit solutions of the boundary 
value problem for a linear system are presented. In Section [3j we present numerical 
schemes for the system of conservation laws (1.1) that converge to the physically 



relevant solution. Second-order schemes are discussed in Section |4j 



2. Theoretical framework 

In the following section, we focus on the so-called boundary Riemann problem, 
which is posed when the Cauchy and Dirichlet data for the mixed hyperbolic- 
parabolic system ( 1.5) are two constant states, U (a;) = U € K m and U;(t) = U; e 
E m . The solution of a general initial-boundary value problem can be build using 
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the solutions of the boundary Riemann problem and standard Ricmann problems 
in the interior, see Goodman [T3] and Sable-Tougeron [21] (Glimm scheme) and 
Amadori [T] (wave front-tracking algorithm). 

2.1. Linear case. 

2.1.1. The solution of the Riemann problem in the linear case. We start by recalling 
the solution of the Ricmann problem obtained by coupling the linear system 



(2.1) 

with the initial datum 
(2.2) 



U t + AU X = U G 



U(0,x) 



u 



x < 
x > 0. 



In (2.1|, A is a constant, strictly hyperbolic m x m matrix, and in (2.2) U + and 



U~ are two given values in R m . 

Denote by Ai, . . . , A m the eigenvalues of A and by R\, . 
right eigenvectors and consider the linear system 



. , R m the corresponding 



(2.3) 



i=l 



which by strict hyperbolicity admits a unique solution (ai, . . . , ct m ). Then the 



solution of (|2.1|)-(|2.2[1 is 

u- 



(2.4) V(t,x) 



U~ 



if x < Ai< 

3 

a.iRi if Xjt < x < \j+it, 
if x > X m t 



, m 



2.1.2. The solution of the boundary Riemann problem in the linear case. We now 
consider the boundary Riemann problem obtained by coupling the linear mixed 
hyperbolic-parabolic system 



(2.5) 



U, e + AVt = eBV' 



Ue 



with the Dirichlet and Cauchy data, 

(2.6) U(t,0)=Ui, U(0,aO = U (aO, 



Vt > 0, x > 0, 



and by taking the limit e — > + . The matrix A in (2.5) is a constant m x m 



matrix satisfying (1.2), and B is another constant, m x m matrix which depends 



on the underlying physical model (we discuss explicit examples later in this paper). 
The data Ui and Uq in (2.6) are constant states in R m . Note that, in general, the 



problem (2. 5), (2. 6) may be ill-posed if the matrix B is not invertible. However, to 



simplify the exposition in the present paper we always choose the data ( 2.6 1 in such 
a way that it is well-posed. 

As mentioned in the introduction, one of the main challenges coming from the 
presence of the boundary is the following: denote by U the limit e — > + of U e , 
then in general the trace of U on the t-axis is not U; , 



(2.7) 



U= lim U(t,x) ^ U;. 

x->0+ 
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More precisely, the relation between U& and U is the following: there is a function 
W : [0,+oo[^ R m satisfying 



(2.8) 



BW = A(U - U) 
W(0) = XJ t , 



lim W(y) 

y— >+oo 



U, 



where we denote by W the first derivative of W(y). A function W satisfying (2.8) 
is called a boundary layer. 



Under the assumption (1.2) and in the case when B is the identity, system (2.8) 



admits a solution W if and only if 



(U; - U) G span(-Ri, ■■■,Rk 



We recall that by (1.2), k is the number of negative eigenvalues of A and as in 
Section 2.1. 1| we denote by R\, . . . , Rk the corresponding eigenvectors. In general, 
when the matrix B is invertible, the system (2.8) admits a solution if and only if 
U ; — U belongs to the stable space of B~ X A (i.e., to the subspace of K m generated 
by the generalized eigenvectors associated to the eigenvalues of B~ X A with strictly 
negative real part). Note that this space depends on the matrix B: this is the reason 
why, even in the simplest possible case (linear system with an invertible viscosity 
matrix), the limit e — > + of (2.5), (2.6) depends on the choice of B. 

In the case when B is not invertible, the analysis in Bianchini and Spinolo [SJ 
Sections 4.2,4.3] guarantees that, in physical cases, if the initial-boundary value 



problem (2.5),( |2.6| is well-posed, then there are k linearly independent vectors 
Ri, . . . ,Rk such that the following two properties hold: first, system (2.8) admits 
a solution if and only if 



(2.9) 



(Uj -U) € span^i, ...,R k 



Second, the vectors Ri, . . . , Rk, Rk+i, ■ ■ ■ , R m constitute a basis of M. m . Specific 
examples with explicit constructions of the vectors Rx, . . . , Rk are discussed later. 
Consider the linear system 



(2.10) 



U, 



k 



i.Ri. 



E 

i=k+l 



atiRi = Uq, 



which by the second property of the vectors Ri , . . . , Rk admits a unique solution 



(«i, . . . ,a m ). The solution U obtained by taking the limit e — >• + of (2.51,(2.6) is 

then 

(2.11) 



if < x < X k t 



U(t,x) = < 



i=l 
k 



Ub + J]] Q4R4 + ctiRi if Xjt < x < \j+\t, j = k + 1, . . . , m — 1 

i=fe+l 



i=l 



I U 



if x > X m t 



where as usual Ai, . . . , A m denote the eigenvalues of the matrix A. Note that this 
construction also works in the case when the matrix B is the identity provided that 
we set 



(2.12) 



Ri—Rj 



Mi = 1, . . . , k. 
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2.2. Explicit computations for the linearized shallow water equations. 

The above constructions are fairly general and abstract. We illustrate them by an 
example, i.e. the linearized shallow water equations of fluid flow (see LeVeque |23j): 



(2.13) 



tk - 

Ut 



uh x 

■gh x 



hu x = 0, 

UU T = 0. 



Here, the height is denoted by h and water velocity by u. The constant g stands 
for the acceleration due to gravity and h, u are the (constant) height and velocity 
states around which the shallow water equations are linearized. 

The physically relevant viscosity mechanism for the shallow water system is the 
eddy viscosity. Adding eddy viscosity to the linearized shallow water system results 
in the following mixed hyperbolic-parabolic system: 



(2.14) 



h t - 



uh x 

gh x 



hu x = 0, 



For the sake of comparison, we add an artificial viscosity to the linearized shallow 
waters by including the Laplacian. The resulting parabolic system is 



(2.15) 



ht- 



uh x 

■gh x 



hu x 



eh r , 



Systems (|2.15|) and (|2.14|) can be written in the form (|2.5|) provided that 
(2.16) A = 



u h 
9 u 



B = B Lap = 



1 
1 



j^EDvisc 



in (2.15) and (2.14), respectively. We will construct explicit solutions for the lin- 



earized shallow water equations (2.13) for the limit of both the eddy viscosity as 



well as the artificial viscosity. For the rest of this section, we specify the parameters 



(2.17) h 
and consider the initial data 



1. 



1. 



(2.18) 



(MXM) = 




and the Dirichlet boundary data 
(2.19) {h,u)(-l,t) 



Ui(t) = (2,1) Vt>0. 



2.2.1. Solution of the Riemann problem. We now apply the construction described 



in Section 2.1.1 to solve the Riemann problem (2.13), (2.18). The eigenvalues of the 



matrix A in (2.16) are Ai = 1 — y/2 < and A2 = 1 + v2 > 0, with corresponding 
eigenvectors 



(2.20) 



Ri 



1 

-V2/2 



R2 



1 

V2/2. 



Hence, the solution of the linear system (2.3 1 in this case is ot\ = 0*2 = 



-1. 
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2.2.2. Solution of the boundary Riemann problem limit of the uniform (Laplacian) 
viscosity. We now apply the construction in Section |2.1.2| to determine the limit 
e — > + of the viscous approximation (2.151 coupled with the Cauchy and Dirichlet 



data (2.181 and (2.19) 



We solve the linear system (2.101 in the case when U; is given by (2.19), Uq 
(3, 1) and (2.12) holds and we obtain a\ — ct^ — 1/2. 



By combining this with the analysis in Section 2.2.1 and by recalling (2.4) 



and (|2.11|), we conclude that the local in time solution of (|2.13|) obtained by taking 



the limit e 
(2.21) 

(h,u)(t,x) 



of (2.151, (2.18), (2.191 is 



(5/2, 1 - V2/4) 
(3,1) 

(2, 1 + y/2/2) 
(1,1) 



if < x + 1< (1 + V2)t 
if x + 1 > (1 + y/2)t and x < (1 
if (1 - \/f)t < x < (1 + y/2)t 
if x> (1 + V2)t 



V2)t 



2.2.3. Solution of the boundary Riemann problem limit of the eddy viscosity. We 
evaluate the limit e — > + of the viscous approximation (2.14) coupled with the 



EDi 



Cauchy and Dirichlet data (|2.18|) and (|2.19|) by applying the construction described 
in Section 
B 



2.1.2 We consider system (2.8) in the case when B is the same matrix 



as in (2.16) and Ui is given by (2.19) and we get 

h- h + 2(u- u) 
— (u — u) 

u(0) 



° 

(2.22) i it 

{ (h,u)(0) = (2,1) 

By imposing the initial datum h(0) 
solution if and only if 

2-h 
l-u 



lim„ 
2, 



(h,u). 



1 one gets that (2.22) admits a 



= span( J Ri) R 1 =( 1,-1/2 



and hence by solving the linear system (2.10) we get in this case ot\ 
and «2 = l/(v / 2 + !)■ By combining this with the analysis in Section 



2.2.1 



2 + 1) 
and 



by recalling \2A\ and (2.11) we conclude that the local in time solution of (2.13) 

obtained by taking the limit e -> 0+ of ( |2.14[ ), ( |2.18[ ), ( |2.19| is 

(2.23) 



(h,u)(t,x) = 




l),(V2 + 2)/(2V2 + 2) 



if < x + 1< (1 + y/2)t 
if a; + 1 > (1 + y/2)t and x < (1 
if (1 - \/2)* < a; < (1 + >/2)t 
if x> (1 + V2)< 



The explicit calculations clearly show that the solution of the linearized shallow 



water equations (2.13) realized as a limit of vanishing eddy viscosity (2.14) differs 



from the solution realized as a limit of the artificial viscosity (2.15) 
the solutions are different near the boundary at x 



In particular, 



whereas they are the same, 
away from the boundary. The height (h) for both solutions is shown in figure [lj 
left. 

2.3. The solution of the boundary Riemann problem for non-linear sys- 
tems. Consider the boundary Riemann problem obtained by coupling the mixed 



hyperbolic-parabolic system (1.5) with the Cauchy and Dirichlet data U e (0,x) = 
U , U e (£, 0) = U;, respectively, where U ,U; 6 E m , and by then taking the limit 
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e — » + . One of the main challenges posed by this problem is establishing the 
relation between the data U; and the trace 

U= Km V(t,x). 

K-J-0+ 



As pointed out by Gisclon and Serre [12] [13], there exists a boundary layer W : [0, 
such that 



(2.24) 



B(W)W = F(W) - F(Z7) 

W(0) = XJ t lim^+oo W(y) = U. 



In the case when the Jacobian matrix Fu satisfies (1.2) and the matrix B is in- 



vertible, the existence of a boundary layer W satisfying (2.24) is equivalent to the 



fact that U; belongs to a suitable stable manifold centered at U. The general case 



when condition ( 1.2 ) is violated (i.e., one eigenvalue of the Jacobian matrix Fu can 
attain the value 0) or the viscosity matrix B is singular is more complicated, but 
it can be treated under suitable assumptions. A characterization of the solution of 
the boundary Riemann problem obtained as limit of e — > + of the viscous approx- 



imation (1.5) is provided in Bianchini and Spinolo [5] under the assumption that 
the Dirichlet and boundary data are sufficiently close. We also refer to Joseph and 
LeFloch [TS] , Ancona and Bianchini [3] and to the references therein for characteri- 
zations of solution of the boundary Riemann problem obtained as limits of different 
viscous approximations. We observe that in general it is difficult to compute the 
solution of the boundary Riemann problem for a non-linear system explicitly. 

3. Numerical schemes. 
3.1. Definition of the schemes. We write down a semi-discrete conservative 



finite difference (finite volume) scheme for the system of conservation laws (1.1) as 
(3-1) |u j -(t) + ^(F J . +1/2 -F J ._ 1/2 )=0. 

Here Uj ss \J(xj) with Xj being the midpoint of the cell [2^-1/2, 2^+1/2]- As 
stated in the introduction, the numerical flux F J+1 /2 = F(Uj,U J+ i) is determined 
by (approximate) solutions of the Riemann problem at the interface Xj + i/ 2 - An 
equivalent expression of the numerical flux (see the book by LeVeque [33]) is 

(3.2) F j+i/2 = - ^ + 2 - — - 2^i+i/2> 

with T> = T>(XJj, Uj+i) being the corresponding numerical diffusion operator. As 
an example, the Roe diffusion operator (see again [23]) is given by 

( 3 - 3 ) V ]+i/2 = 1 A i+1 / 2 1 i2T +1/2 [U] i+1/2 , 

where A and R are matrix of eigenvalues and eigenvectors of the Jacobian Fu 
evaluated at a suitable average state. Also, here and in the following we use the 
notations 

T 

Clearly, the above numerical diffusion operator does not incorporate any informa- 
tion about the underlying viscous approximation. Hence, the approximate solutions 
generated by schemes such as the Roe scheme may not converge to the physically 
relevant solution, given as a limit of the underlying viscous approximation as the 
viscosity parameter goes to zero. An illustration is provided in figure[T] right. Here, 



a j+i/2 - 7) ■ Hj+1/2 - a J+i ~ a r 



10 



SIDDHARTHA MISHRA AND LAURA V. SPINOLO 



we present results obtained by approximating the linearized shallow water system 



(2.131 with the Roe (Godunov) scheme. The exact solution, computed in (2.23) as 



the limit of vanishing eddy viscosity (2.141 is also shown on the left. As shown in 



the figure, the Roe scheme does not converge to this physically relevant solution as 



the numerical viscosity (3.3) is very different from the eddy viscosity in (2.14) 



In this paper, we consider a different paradigm for numerically approximating the 
initial-boundary value problem. The main difference from the standard schemes lies 



in the choice of the numerical flux F in (3.1 ). It combines the following ingredients 



3.1.1. Entropy conservative fluxes: We assume that ( 1.1 ) admits an entropy-entropy 
flux pair (S, Q) and by following Tadmor 28J we define an entropy conservative flux 
for (1.1) as 



Definition 3.1. A numerical flux F* +] y 2 = F*(Uj, Uj+i) is defined to be entropy 
conservative for entropy S if it satisfies 



(3.4) 



m- +1/2 F* +1/2 = [*] i+1/2 



for every j . Here, V = Su is the vector of entropy variables and "J" = V T F 
the entropy potential for the entropy function S and entropy flux Q. 



Q is 



The existence of entropy conservative fluxes for system of conservation laws is 
shown in Tadmor [28 and explicit examples of entropy conservative fluxes are 
summarized in Fjordholm, Mishra and Tadmor |11) . 

It is known (see Dafermos [5J Chapter 7]) that if (S,Q) is an entropy-entropy 
flux pair (S, Q) with S strictly convex, then F;/Uy is symmetric and Uy is sym- 
metric and positive definite. In the following, we also assume that the entropy is 
dissipative, namely that 



(3.5) 



(BU v £,0>0 V£e 



This condition is satisfied in physical cases. In particular, it is satisfied in all the 
cases we discuss in the following. 



3.1.2. Numerical diffusion operator: Let (1.5) be the underlying mixed hyperbolic- 
parabolic regularization of the hyperbolic equation ( |1.1[ ). We choose a numerical 
diffusion operator, 



(3.6) 



■D* /2 := V{U jt V j+1 ) = c max B(U j+1/2 )[U] 



3+1/2- 



Here, B is the viscosity matrix in the parabolic regularization (1.5) evaluated at 
some suitable averaged state C/ J+1 / 2 and 

(3.7) c max (i)=max|A™ ax |, 

3 J 

with A™ ax being the largest eigenvalue of the Jacobian Fu at a given state XJj . 
3.1.3. Correct numerical diffusion (CND) scheme. We choose the numerical flux 

(3.8) I 



' 3+1/2 



T71* 1 T-) + 

*i+l/2 - 2 J+l/2- 



Here, F* +1 y 2 = F* (Uj,Uj + i) is an entropy conservative flux (3.4) for the system 



(1.1) and the numerical diffusion operator T> * +1 y 2 i s defined in (3.6). The semi- 



discrete scheme (3.1) with numerical flux (3.8) has the following properties 



ACCURATE NUMERICAL SCHEMES 



11 



Theorem 3.1. Assume that the system (1.1) is equipped with the entropy-entropy 
flux pair (S,Q) which is dissipative in the sense of (3.5). Then, the scheme (3.1) 
with numerical flux (3.8) satisfies 

(i.) a (local) discrete entropy inequality (discrete version of the entropy inequal- 
ity ( 1.4 \) of the form 

(3.9) Jt S{V ' ){t) + h (^ +1/2 " ^ +1/2 ) " °' 

with a numerical entropy flux Q that is consistent with the entropy flux Q. 
(ii.) The scheme is first-order accurate and the equivalent equation is 

(3.10) Uf x + F(U Aa % = (£(U A *)U Aa % + 0(Ax 2 ). 

Proof. Multiplying both sides of the scheme ( |3.1[ ) by the entropy variable Vj , we 
obtain 



(3.11) 



2Aa 



[V]J +1/2 F* +1/2 + [V]J_ 1/2 F*_ 1/2 



Ti 

- ~ ([V]J +1/2 S(U(V, +1/2 ))U v (V i+1/2 )[V] i+1/2 
(mJ-i/ a BCU(V i _ 1/2 ))U v (V j _ 1/a )[Vl i -i / 



4/^ V 11 Jlj-1/2"V^ V" j-l/2) J" V V » 3-1/2/1 * II j — 1/2 

Here, we have introduced the numerical flux 

Q j+ i/2 ■= vJf* +1/2 - (vJ +1/2 s(u(v J+1/2 ))[u] j+1/2 ) , 

with Qj-i/2 defined analogously. Also, we have introduced the average V J+1 / 2 
satisfying 

(3-12) [U] i+1/2 = U v (V i+1/2 )[V] i+1/2 . 



By using the definition of the entropy conservative flux (3.4), the term T\ in (3.11) 
can be simplified as 

Tl = 2zL (l V lJ+i/ 2 F ^i/2 + [V]J_ 1/2 F*_ 1/2 

= ^(*i+i/a-*;-i/2) 
Substituting the above expression of T\ in ( |3.11 ) yields 
d , 1 



-5(U j) + ^^ + i /2 -0^i /2 
( 3 - 13 ) " ([V]J +1/2 S(U(V, +1/2 ))U v (V i+1/2 )[V] i+1/2 ) 

- ([V]J_ 1/2 B(U(V,_ 1/2 ))Uv(V,_ 1/2 )[V] j _ 1/2 



with 

Q,+i/2 := vJ +1/2 F* +1/2 - * j+1/2 - (vJ +1/2 B(U(V j+1/2 ))[U] J . +1/2y 
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and Qj-1/2 defined analogously. Note that Q{a 1 a) = Q(a) from the definition 
of the entropy potential 'J. The discrete entropy inequality (3.9| follows from 



assumption (3.5 ) 



The equivalent equation (3.101 is a simple consequence of Taylor expansion and 
reveals the first order accuracy of the scheme. □ 

Thus, the proposed numerical scheme is entropy stable under reasonable hy- 



potheses on the system (1.1). Furthermore, the equivalent equation (3.10) shows 



that the numerical viscosity of this scheme matches the underlying physical viscos- 
ity operator in (|1.5[ ) at leading order. Hence, we claim that the scheme (3.1) with 
numerical flux (3.8) incorporates the correct numerical dissipation and term it as 



the CND scheme. 



The Dirichlet boundary conditions for (1.1) are imposed weakly by setting 



(3.14) U (t) = U,(t). 

This amounts to setting the Dirichlet data as the value in the ghost cell [x—i/a, ^1/2] ■ 
The semi-discrete scheme (3.1) is integrated in time with the SSP-RK2 time 
integrator: 



(3.15) 



ir = u' ; 

IT* = U* 

1 



- At£(U"), 



U n+1 = -(IT 
2 V 



that approximates the ODE system 

d 
dt 



(3.16) 



IT*), 
U(i)=£(U(t)), 



for the unknowns U = {Uj},-, defined by the scheme (3.1 ). 



3.2. Linear systems: We illustrate the finite difference scheme (3.1) for a linear 



system, i.e for (1.1) (and the parabolic regularization (1.5)) with 
(3.17) F(U)=^U, E(U)=B, 

for A,B given (m x m)-matrices. As pointed out before, if S(XJ) = ^XJ T SXJ is a 
strictly convex entropy, then the matrix SA is symmetric. Following Fjordholm, 
Mishra and Tadmor |11|). we define the corresponding entropy conservative flux as 

(3.18) 

We consider the following specific example: 



F J + l/2 - A\J J + 1 /2- 



3.2.1. Linearized shallow water system. The linearized shallow water system (2.13) 



is considered. We assign the data (2.18) and (2.19). The computational domain is 



[—1,1] and we use open (Neumann type) boundary conditions at the right boundary 
x = 1. 



CND scheme (3. 



The numerical solutions computed with the standard Roe scheme (3.1) and the 
are shown in figure [2] As we are interested in computing the 



physically relevant solutions of the linearized shallow water equations, obtained as 
a limit of the eddy viscosity (2.14), we also plot the exact solution computed in 
(2.23) for comparison. Both the numerical solutions are computed with a 1000 
mesh points. 
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The results in figure [2] clearly show that the Roe scheme does not converge to 
the physically relevant solution (2.23|. On the other hand, the solutions computed 



with the CND scheme approximate the physically relevant solution (2.23) quite 
well. There are some small amplitude oscillations in the height with the CND 
scheme. This is a consequence of the singularity of the viscosity matrix B in this 
case. The experiment clearly shows that incorporating explicit information about 
the underlying viscous mechanism in the numerical diffusion operator results in the 
approximation of the correct solution. 



1 



-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



(a) Height (h) 





-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



(b) Velocity (it) 



Figure 2. Solutions of the linearized shallow water equations 



(2.13) with initial data (2.18) and boundary data (2.19) computed 



with the Roe and CND schemes with 1000 mesh points. The exact 



solution computed in (2.23) is provided for comparison. 



3.3. Nonlinear Euler equations. In one space dimension, the Euler equations 
of gas dynamics are 

Pt + (pu) x = 0, 

(3.19) (pu) t + (pu 2 +p) x = 0, 

E t + ({E+p)u) x = 0. 

Here, p is the fluid density and u is the velocity. The total energy E and the 
pressure p are related by the ideal gas equation of state: 

(3.20) E = + \pu\ 

7 — 1 2 



Here, c = ^/-y is the sound speed 
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with 7 > 1 being a constant specific of the gas. 
The system is hyperbolic with eigenvalues 

(3.21) Ai = u — c, A2 = u, A3 = u + c. 

jV_ 

p 

The equations are augmented with the entropy inequality 

(^), + (^o. 

with thermodynamic entropy 

s = log(» -7log(». 

The compressible Euler equations are derived by ignoring kinematic viscosity 
and heat conduction. Taking these small scale effects into account results in the 
compressible Navier- Stokes equations: 

Pt + {pu)x = 0, 

^ 2 ^ (PU)t + {PU 2 + P) X = VU; 

E t + ((E+p)u) x = v 
Here, 9 is the temperature given by 

V 



XX ', 



2 



(7-l)P 

The viscosity coefhcient is denoted by v and k is the coefficient of heat conduction. 

For the sake of comparison, we add an uniform (Laplacian) diffusion to obtain 
the compressible Euler equations with artificial viscosity: 

Pt + (pu) x = ep xx , 

(3.24) (pu) t + {pu 2 + p) x = e(pu) xx , 

Et + ((E + p)u) x = eE xx . 



To evaluate the limit solution of ( 3.23[ ), we construct a numerical approximation 



by discretizing the mixed hyperbolic-parabolic systems (3.23) and ( 3.24[ ) for a fixed 



and very small value of the viscosity coefficient. We do so by the (semi-discrete) 
finite difference scheme 

(3-25) l U i(*) + i( F W- F ^/ 2 )=^ 



Here, the numerical flux is the entropy conservative flux (3.4) (see Ismail and 
Roe PS]) 

TT* — fTT 1: * TP 2 '* "IT 3 '* 1~T 

r j + l/2 ^ r i + l/2' r j + l/2' r J+1/2J ' 

(3.26) i +1 / 2 ~ Wi+WWi+V* - (zij^ + F ^V2' 

3,* _ 1 CgDj+vg ( 7 + 1 ( z 3)j + i /2 a,. \ 

*^ 2 2(zr), +1/2 ^7-i(^ + i/2 J+1/ 7 
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with parameter vectors 
(3.27) (21,22,23) 

' " p y p 

The logarithmic mean of any quantity a, defined on the mesh, is denoted by 




«+l/2 = 



a j+ i 



>i +1 ' 2 log(a J+1 )-log(a i )' 
We define the numerical diffusion operators by setting 



2> J =[2>),2> i 2 ,Zf] T 



(3.28) 



and 



(3.29) 



V 2 = v{u 



2uj + M7-1) 



2) i = 2^ +1 " 2 " 



2 1 u 2 



K(9 j+1 - 29 j + 9 j _ 1 ). 



V j =U j+1 -2U j + V 



i-i) 



for the compressible Navier-Stokes equations (3.23) and the Euler equations with 



artificial viscosity (3.24), respectively. 



As an example, we consider both (3.24) and (3.23) in the domain [—1,1] with 
initial data 



(3.30) 



{p ,u ,p ) 



(3.0,1.0,3.0), ifx<0, 
(1.0,1.0,1.0), ifx>0. 



We impose open boundary conditions at the right boundary and Dirichlet boundary 
conditions at the left boundary with boundary data 



(3.31) 



(p(-l,t),u(-l,t),p(-M)) = (2.0,1.0,2.0). 



We set v = k = e. The results for the finite difference scheme approximating the 



uniform viscosity (3.24) and the physical viscosity (3.23) are presented in Figure 



[3] The figure shows that the there is a clear difference in the limit solutions of this 
problem, obtained from the compressible Navier-Stokes equations (3.23) and the 



Euler equations with artificial viscosity (3.24). The difference is more pronounced in 



the density variable near the left boundary. Both the limit solutions were computed 
by setting e = 10 -5 and on a very fine mesh of 32000 points. 

Remark 3.1. The above example also illustrates the limitations of using a mixed 



hyperbolic-parabolic system like the compressible Navier-Stokes equations (3.23). 
In order to resolve the viscous scales, we need to choose Ax = O {-), with e being 
the viscosity parameter. As e is very small in practice, the computational effort 
involved is prohibitively expensive. In the above example, we needed 32000 points to 
handle e = 10~ 5 . Such ultra fine grids are not feasible, particularly in several space 
dimensions. 
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(a) Density (p) 
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-PHYVISC 
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(b) Velocity (u) 



j 



■LAP 

-PHYVISC 



-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



(c) Pressure (p) 



Figure 3. Limit solutions of the compressible Euler equations 



(3.19) with initial data (3.301 and boundary data (3.31). The limits 



of the physical viscosity i.e compressible Navier-Stokes equations 



(3.23) and the artificial (Laplacian) viscosity (3.24) are compared 



3.3.1. CND scheme for the Euler equations. The CND scheme (3.1) for the Euler 
equations (3.19) is specified as follows: the entropy conservative flux in (3.8) is 
given by (3.26) and the numerical diffusion operator in (3.8) matches the kinematic 
viscosity and heat conduction of the compressible Navier-Stokes equations since it 
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is defined by setting 



^ , j+l/2) :D i+l/2' I> i+l 



(3.32) 



V j + l/2 



^1+1/2 = (^max ( //, + 



/2 



n 3 

^3 + 1/2 



max kt 7 - + 



1 - 



We discretize the initial-boundary value problem for the compressible Euler equa- 



tions (3.19) on the computational domain [—1, 1] with initial data (3.301 and Dirich- 
let data ( 3.31[ ). The results with the CND scheme and a standard Roe scheme are 
shown in figure [4] We present approximate solutions, computed on a mesh of 1000 
points, for both schemes. Both the Roe and the CND schemes have converged at 
this resolution. As we are interesting in approximating the physically relevant solu- 
tions of the Euler equations, realized as a limit of the Navier-Stokes equations, we 
plot a reference solution computed on a mesh of 32000 points of the compressible 
Navier-Stokes equations (3.231 with k — v = 10~ 5 . The figure shows that the Roe 



scheme clearly converges to an incorrect solution near the left boundary. This lack 
of convergence is most pronounced in the density variable. Similar results were 
also obtained with the standard Rusanov, HLL and HLLC solvers (see the book by 
LeVeaue[2"3"] for a detailed description of these solvers). 

On the other hand, the CND scheme converges to the physically relevant solu- 
tion. There are slight oscillations with the CND scheme as the numerical diffusion 
operator is singular. However, these oscillations do not impact on the convergence 
properties of this scheme. Furthermore, the CND scheme is slightly more accurate 
than the Roe scheme when both of them converge to the same solution (see near 
the interior contact). 



4. Second-order CND schemes 

The CND scheme, described in the last section, was first-order accurate in space. 
Consequently, it approximated shocks and contact discontinuities with excessive 
smearing, particularly on coarse meshes. We can improve the resolution of numer- 
ical schemes by constructing second-order accurate schemes. 

To this end, we reconstruct the cell averages Uj of the unknown to a piecewise 
linear function given by 

U' 

(4.1) p j (x):=U j + -^(x-x j ). 

The numerical derivative U' is chosen to be non-oscillatory by limiting the slope, 
i.e. setting 

(4.2) Uj- = minmod(U i+ i - Vj,Vj - Uj-i), 
with the minmod function defined as 

sgn(a) min{|a|, if sgn(a) — sgn(b), 



(4.3) minmod(a, b) 



0, otherwise. 
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Figure 4. Approximate solutions of the compressible Euler equa- 



tions (3.19) with initial data (3.30) and boundary data (3.31). We 



compare the Roe and CND schemes on 1000 mesh points with 
a reference solution of the compressible Navier-Stokes equations 



(3.23) with k = v= 10" 



Other limiters like the MC and Superbee limiters can also be chosen (see the book 
by LeVeque[23j for the corresponding definitions). We need the cell interface values 



(4.4) 



PiO. 



u: 



Pj(%-l/2)- 



With these reconstructed values, we modify the numerical flux ( |3.8[ ) by setting 
(4.5) F J+1/2 



F j+i/2 - 2^J+l/2) 
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with 



(4.6) V. ]+1/2 = P(U+ UJ +1 ) = c m£ 

where the constant c max is the same as in 



B(U J+1/3 )(U7 +1 -Ut), 
3.7|). Note that the on/y difference be- 



tween the flux (3.8 1 and the flux (4.5 1 lies in replacing the difference in cell averages 



in the numerical diffusion operator in (3.6) with the difference in the corresponding 



reconstructed edge values in (4.6). The overall scheme (3.1) with numerical flux 



(4.5) is (formally) second-order accurate as the entropy conservative flux F* is sec- 
see Tadmor 



ond order accurate (sec Tadmor 28 ) and the difference in the numerical diffusion 
operator is a difference of second-order reconstructed values, see Fjordholm, Mishra 
and Tadmor [11] for a proof of the order of accuracy of schemes constructed with 
numerical fluxes like (4.5). 



We test this second-order scheme (3.1), (4.5) for the compressible Euler equa- 



tions. Let the computational domain be [—1, 1] with initial data (3.30 ) and Dirichlct 



data (3.311 



The scheme (3.1 1 is specified as follows: the entropy conservative flux in numer- 
ical flux (4.5) is given by (3.26). The numerical diffusion is 



V 



J+l/2 



^J + l/2' 



T) 2 T) 3 

^'+1/2)^+1/2 



(4.7) 



1/2 



D 3 



1/2 



max 



max 




«) 2 ) 



with u ,0 being obtained from the reconstructed conservative variables. The 



overall scheme (integrated in time with the SSP RK2 time stepping (3.15)) is termed 
as the CND 2 scheme. 



We compute approximate solutions of the Euler equations with initial data (3.30 ) 



and boundary data (3.31 ) using the CND and CND2 schemes and show the results, 



obtained on a mesh of 200 points, in figure [5j The result shows that both the 
first and second order CND schemes approximate the physically relevant solution, 
computed as the limit of the compressible Navier-Stokes equations, quite well. The 
first-order scheme smears the discontinuities as well as generates oscillations. On 
the other hand, the second-order scheme is clearly sharper at discontinuities. Fur- 
thermore, it reduces the oscillations considerably. 



5. Conclusion 

We consider the initial-boundary value problem for systems of conservation laws 



(1.1). Since the work by Gisclon and Serre [T2j [13] it is known that the solutions 
of the initial boundary value problem depend on the underlying viscous approxi- 



mation (1.5). Different choices of viscosity operators can lead to different solutions 



for the limit system of conservation laws (1.1). These results hold for both linear 



as well as non-linear systems. Even 2x2, strictly hyperbolic, symmetrizable linear 



systems like the linearized shallow water equations (2.13) show this behavior. 



This dependence of solutions on underlying small scale effects suggests that one 
should discretize the viscous approximation (1.5) directly. However, this is very 



expensive computationally on account of very low values of the viscosity parameter. 
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(c) Pressure (p) 



Figure 5. Approximate solutions of the compressible Euler equa- 
tions (3.19) with initial data (3.30) and boundary data (3.31). We 



compare the CND and CND2 schemes on 200 mesh points with 
a reference solution of the compressible Navier-Stokes equations 



(3.23) with k = v = 10~ 5 . 



Therefore, we need to design numerical schemes for the system of conservation laws 



(1.1) that converge to the physically relevant solutions i.e the limit of solutions 
of (1.5) as e — > 0. Unfortunately, existing numerical schemes like the standard 



Godunov, Roe and HLL schemes might converge of the physically incorrect solution 
of the initial-boundary value problem. 



In this paper, we design a conservative finite difference scheme (3.1) with a 
based on the following two ingredients: 



numerical flux (3. 



entropy conservative fluxes (3.4) 
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numerical diffusion operators (3.6) 



Information about the underlying viscous approximation (1.5) is explicitly incor- 
porated into the choice of the numerical diffusion operator. The resulting entropy 
stable schemes are shown (numerically) to converge to the limit solution, obtained 
from the underlying viscous approximation. Thus, we provide a numerical frame- 
work for computing solutions of the system of conservation laws that require explicit 
information about the underlying small scale effects. To the best of our knowledge, 
this is the first time such schemes have been constructed in the context of initial- 
boundary value problems. 

We present a set of numerical experiments for both the linearized shallow wa- 
ter and nonlinear Euler equations to demonstrate that our numerical schemes do 
converge to the limit solutions of the underlying eddy viscosity or Navier-Stokes 
viscosity, respectively. Second-order schemes are constructed and are shown to be 
superior to first-order schemes in terms of accuracy as well as in suppressing oscil- 
lations that might result from a lack of viscosity in some conservative variables. At 
the same time, these second-order schemes also converge to the physically relevant 
solutions. 

We concentrated on Dirichlet boundary conditions in one space dimension in this 
paper. Extensions to several space dimensions and to other interesting boundary 
conditions will be considered in a forthcoming paper. 
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